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. Abstract 

We present an algorithm to reconstruct gray scale images corrupted 
by noise. We use a Bayesian approach. The unknown original image 
is assumed to be a realization of a Markov random field on a finite 
00 ! two dimensional region A C Z 2 . This image is degraded by some 

noise, which is assumed to act independently in each site of A and to 

■ have the same distribution on all sites. For the estimator we use the 
mode of the posterior distribution: the so called maximum a posteriori 
(MAP) estimator. The algorithm, that can be used for both gray-scale 
and multicolor images, uses the binary decomposition of the intensity 
of each color and recovers each level of this decomposition using the 
identification of the problem of finding the two color MAP estimator 
with the min-cut max-flow problem in a binary graph, discovered by 

. Greig, Porteous and Seheult (1989). Experimental results and a 

detailed example are given in the text. We also provide a 

■ web page where additional information and examples can be 
found. 
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1 Introduction 



We consider here the problem of image reconstruction. Suppose a multicolor 
or gray-scale picture is subjected to noise and an observer has access only to 
this corrupted version. How can she estimate the original picture? 

The analysis of this kind of problem has attracted a lot of interest and 
many approaches have been considered ([G]). 

One of the methods proposed is the so called MAP estimator. In this 
method one assumes that the original image is a random realization of a 
Markov random field that has been corrupted by some site independent noise. 
One assumes that the distribution of the field (a priori distribution) is known 
as well as the distribution of the noise, that is, the conditional distribution of 
the observed image given the original one. The MAP estimator is the image 
that has the largest probability of have produced the observed one. This is 
the mode of the posterior distribution. 

In multicolor images the literature proposes to use the so called Potts 
model as a priori distribution. Roughly speaking, this model is a measure 
on the set of images that gives more weight to images that have neighboring 
pixels of the same color. This choice has a number of advantages and dis- 
advantages. An important disadvantage is the fact that the algorithms used 
to compute the MAP estimator operate in exponential time in the number 
of pixels ([GJ]). Our main point in this paper is to propose an alternative 
a priori measure with the property that the computation time of the exact 
MAP estimator is polynomial in the number o pixels. Using this approach 
we have produced an algorithm and a program that reconstruct dirty images 
in polynomial time in the number of pixels. 

In the remaining of this section we explain these ideas in some detail. 
In the next section we present some experimental results obtained from the 
implementation of our method to some images. The next section is more 
technical in nature. It explains why our MAP estimator can be computed ef- 
ficiently (polynomial time) and provides a proof for the theoretical result pre- 
sented below. The question of how to obtain estimators for some important 
parameters associated to the observed image is discussed in the appendix. 
We close this paper with some final remarks. 

In order to motivate the discussion we have to introduce some notation. 
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Assume the image is a point of {1, • • • , c} A , where A is a finite subset of Z 2 
(e.g. a square) with N sites. The image associates to each site or pixel in A 
one of c possible colors. Typically values for c are c = 2 8 or c = 2 24 . Assume 
1<K = log 2 c E Z and let X = (^) i=1 ,... iiV for X, = (X/, X?, . . . , X? ) e 
{0, V « G A, X > 1, represent the true unknown (random) image. 
Denote the observed image by Y_ = (Yi) i=1: ___ :N . We denote the space of 
pictures {{0, 1} X } A by S^. In this notation we have c = 2 K possible colors 
in each pixel. Each Xi may correspond to a single binary number which 
gives the intensity of black in pixel i as J2k=i Xf2~ k (in the case of gray- 
scale picture) or it may correspond to three binary numbers, each giving the 
intensity of one of the three basic colors. In this second case, for instance, we 
could have K = 24 and Y^=i+c8 Af 2 _fc , for c = 0, 1 and 2 giving the intensity 
of, respectively, Red, Green and Blue at site %. To simplify the discussion we 
assume, without loss, the first case since our approach in the second one is to 
reconstruct each color separately in order to reconstruct the whole picture. 

In the Bayesian setting we assume (a) that the original picture A is 
random with a known distribution which is called the a priori measure and 
(b) that we know how to model the noise. 

To start we consider the noise. We fix an e > 0. Conditional on A, at 
each bit k of each pixel i, independently, the observed value in bite k of pixel 
i, called Yf, is equal to the true value X\ with probability 1 — e or, with 
probability e, it corresponds to the switched value: 




(1) 



Hence 



P(Y\2Q 



n n 



(2) 



ieAke{l,...,K} 



1 
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h = log 



(1-*) 
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and Z h is the normalization constant. 
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Note that this hypothesis is not the same as the one in [FFG] where each 
pixel was either observed correctly with probability 1 — e or chosen uniformly 
among the other c — 1 colors. 

Now we define the MAP (maximum a posteriori) estimator. 

Denote by \i the a priori measure, that is, the distribution of X_. Then 
the MAP estimator after observing the image Y , X — X(Y) G T^k is any 
image which maximizes the posterior distribution P(X\Y) oc P(Y \X)fi(X), 
that is 



P{2L\Y) = maxP(X\Y) (5) 

Now we discuss the a priori measure. We want to consider here the 
situation on which very little is known about the original picture before 
the observation besides the information that it is some kind of real photo, 
perhaps taken by a satellite, as opposed to being a photo of something like 
a geometrical drawing, a cubist oil painting or a cell of a cartoon picture. 
We suppose that the only prior information available is that the original 
picture is locally smooth. By this we mean that the measure should be such 
that neighbor pixels on the picture are more likely to have colors which are 
near in some sense. More precisely we assume that there exists a real valued 
function H which indicates how smooth is an image and take the a priori 
measure to be: 



P -PH(X) 

M20 = — ^- (6) 

where Z = J2x e~P H ^—\ with the sum taken over all possible images, is 
a normalization constant and (3 > is a real parameter. This parameter 
measures the tendency to be smooth since if (3 is large [i is concentrated on 
images with small values of H while if f3 is small /i is close to the uniform 
distribution on S^, the set of all images. The motivation for this formula, 
in particular the minus sign in front of /3, comes from statistical mechanics 
where it would be called Gibbs measure, j3 would be the inverse temperature 
and H would be the Hamiltonian (or Energy) function. 

Plugging ([|) into (|5]), and taking logarithms, we get that the images which 
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maximize the posterior distribution are those which maximize 

(3H(X) + hY: E l{*?=i?>£) ( ? ) 

The image that maximizes this expression makes the best compromise be- 
tween being globally smooth, regulated by /3H, and agreeing as much as possi- 
ble with the observed image which is regulated by hY^ieA J2ke{i,...,K} l{x fe =y. fc }(H)- 
Of course there exists only one relevant parameter in the maximization prob- 
lem, say h//3. 

Although any strictly positive measure \x on can be written as above 
for some H, and thus @ can be safely assumed for a generic a priori 
measure, we are interested only in the case on which H is both simple and 
reasonable as a measurement of smoothness. This simplicity requirement 
will basically mean the assumption that H has only local dependence as 
follows 

H(X)= E d{X h X j ) (8) 

<i,j> 

where d(Xi,Xj) is a measure for the distance between the colors at pixel % 
and j and the sum is taken over all pairs of neighbor pixels in the picture. 
That is, < i,j >= {(i,j) : \i - j\ = 1}. 

The choice of d(Xi,Xj) depends on what kind of local properties one 
expects in the original picture. 

One choice which is common in the literature (see, for instance [FFG] 
and its quotations) is 

H P (X) = J2 (9) 

<i,j> 

where the sum is taken over all pairs of neighbor sites in the lattice. 

This choice is related to the Potts model in statistical mechanics (see for 
instance [M]) and explains the subscript. It is a model with very interesting 
properties that corresponds, in the two-color case, to the Ising model [MW]. 
To find a solution to @ we have to find what is called in the physics literature 
a ground state for the Ising model with random magnetic field. This field is 
induced by the observed image. 
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Most of the interesting properties of those statistical mechanics models, 
like phase transition, appear in the so called thermodynamic limit, the limit 
on which the size of the system grows to the whole lattice ([R], [MW]). 
Even though in our discussion the lattice size is kept fixed we will need in 
section 4 some exact results on the thermodynamic limit in the two-color case 
(Ising model). Also some finite size considerations like the effect of boundary 
conditions (in our case we chose free boundary conditions) may be important. 

The measure defined by (|B|), with (^) plugged in, would describe images 
on which neighbor pixels tend to be equal but if two pixels have different 
colors the cost of this interface does not depend on how different they are. 
If two pixels have different colors the most likely is that each one belongs to 
a single-color region separated by a sharp interface. 

We could also try to represent the situation on which this is not neces- 
sarily the case choosing a finer notion of distance between X% and Xj. Two 
somewhat natural choices would be 



Hi(2Q 



E 

<i,j> 



K 



k=l 



2 k 



(10) 



or 



H2(K) = - £ (E(xf-x^)l) • (ii) 

<i,j> \fc=l Z / 

Note that, for any choice of H, the MAP problem is well posed since 
we could, at least in principle, check the finitely many values in the right 
hand side of (|?|) and choose a image that maximizes it. But since |E#| = 
2 K with something like K = 24 and N = 400 x 600 this approach is not 
computationally feasible. The problem with all the above choices is that 
it is not known how to do better than this time-consuming maximization 
by inspection in the multicolor case (K > 1) and thus the problem is, in 
practice, not solvable. 

We propose another choice for H, intermediate between the Potts and the 
other two mentioned above, which respects nicely the notion of smoothness 
which led to the choice flIU| ) but nevertheless induces to a polynomial time 
maximization problem. This choice is: 
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h{x) = EEI^-4lrEE V^xn ^ (12) 

k=l <i,j> k=l <i,j> 

With this, the problem of finding the 2 x -color image which maximizes 
the posterior distribution P(X\Y) oc P(Y\X)ii(X) is decomposed into K 
binary color maximization problems. Namely one has to solve 

P(X k \Y k ) = maxP(n|20M2O (13) 

for each k, 1 < k < K where Y_ k = {Y^}^^} is the k-th component of the 
observed image and fi k is the Gibbs measure defined on the space of binary 
color images, S 2 , by @ with 

H\X k ) = Y, hx^xf} L (14) 

<i,j> 

The K color image given by F is a solution of Each one of these K 
binary images is called a layer. The weighted sum of the MAP for each layer 
gives our (gray-scale) image estimator: 

k=l 

On the multicolor case one solves a gray-scale problem for each basic 
color. 

Each binary problem can be solved in polynomial time using the results 
by Greig, Porteous and Seheult who reformulated it as one involving find- 
ing a minimum cut on a capacitated network [FF] for which there exist fast 
algorithms. These ideas are presented briefly in the next section. From the 
statistical point of view the approach is different whether the parameters (3 
and e are known or not. A truly Bayesian approach would associate an a 
priori measure to each one of the parameters f3 and e. We leave this alter- 
native to future work. Another possibility is to estimate these parameters 
from the observed image using classical frequentist analysis. This is pos- 
sible under the hypothesis on the noise and on the original image being a 
sample of a Gibbs measure for the Potts Hamiltonian (^]) (which is always 
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the case after the decomposition in binary colors) using exact results for the 
two-dimensional Ising model obtained by Frigessi and Piccioni [FP]. We have 
produced an algorithm based on [FP] for estimating the parameters e and (3. 
This is explained in section 4 below. 

Once one has a fast algorithm to reconstruct images it is natural to ask 
what is the effect of iterating the whole process. More precisely, what hap- 
pens if one takes the reconstructed image and applies the method again, 
perhaps updating the value of hj '/3? 

Note that to apply the method again for the (already) reconstructed 
image is equivalent to assume that the reconstructed image could be thought 
as obtained from some original picture which was chosen with respect to fi 
and then subjected to noise. Even though it is not difficult to verify that this 
assumption is false the question is interesting both from the mathematical 
and from the applied point of view. On the mathematical side one has a 
mapping, Y_ i— > X(Y), from into itself, a discrete time dynamical system, 
and it is natural to ask about iteration properties. On the practical minded 
side one can ask about the effect of iteration on the quality of reconstruction 
even if this can only be judged subjectively. 

If we decide to update hj f3 before each iteration the natural thing to do 
would be to chose a larger value at each step since we assume the procedure 
did a good job and removed some of the noise, therefore X(Y) would be more 
reliable than Y_. It would then be natural to increase h, without changing /3. 
A natural procedure to find each updated value is to use again the estimators 
given by Frigessi and Piccioni ([FP]). In doing this we find, in fact, that this 
parameter increases, after one iteration, as expected. 

What is rather surprising is that once one tries this iteration procedure, 
either keeping h/(3 fixed or increasing it, at each step, one finds that it has 
no effect et all. The twice reconstructed image is exactly equal to the once 
reconstructed one. In other words, X as a function from into itself has a 
fixed point in each once- reconstructed image. 

More precisely write F(h/p)(X) = X_ for the function from into itself 
defined in (|^). Then we have the following 

Proposition: F S (X) = F t (F s (X)) for all t > s. 

We prove this proposition in the next section after reviewing some results 
on networks and the connection with the maximization problem considered 
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in this note. 

Given this result a natural question arises. Is this property not true in 
general! Suppose one has a random function assuming values in some space 
S with some unknown parameter s which itself belongs to S chosen according 
to some a priori measure (in our case, s corresponds to the original image). 
Define the MAP estimator as usual to obtain a function from S into itself. 
Does it always satisfy the this fixed point property? One could argue that 
since all the information about the original image which is contained in the 
observed one should be still present in the reconstructed image then this 
iteration should give no further results and the twice reconstructed image 
should always be equal to the once-reconstructed one. As it turns out this is 
not the case and it is not difficult to find counterexamples. 

2 Experimental Results 

We developed a program (map) in C to reconstruct images according to the 
method proposed. This program reads an image in the portable bit map 
format for true-color (.ppm) or gray-scale (.pgm) pictures. The intensity of 
each color in each pixel (or the single gray scale intensity there) is written as a 
(eight bits) binary number. The program then constructs a graph for each one 
of those layers (as defined after equation ([14])), construct the corresponding 
graphs and finds the minimum cut for each one of them using a standard 
implementation of the Ford Fulkerson algorithm. The program uses integer 
arithmetic scaled by a factor of 10000. The values for (3 and H are supplied 
by command line switches. It is also possible to supply a bit mask to select 
which planes to reconstruct. The program is portable and it has been tested 
on SunOS, Solaris and Linux. 

The main problem in developing the program was the large amount of 
memory to hold the image, the graph and the auxiliary data structures. The 
solution was to use a compact representation of the graph, namely a matrix 
where each cell has a pixel and the values for the flow in each of 5 directions 
(4 neighbors plus s or t). This avoided the explicit representation of the 
edges and the extra space for the image itself. 

A series of tests were made on a 700Mhz Pentium III machine, with 
128Mb of main memory, running Linux Debian potato 2.2, kernel version 
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2.4. The running times for restoration of six sample images are presented in 
table |l|. These images are shown in tables |2|, |3| and f|. 

The images were obtained from a picture shot with a digital camera 
and converted to 8 bit gray-scale. The final result was saved in the file 
"f ish.pgm". Two aditional files were obtained by clipping the image to a 
rectangle with half of the area ("f ish2.pgm") and to another one with a 
quarter of the total area ("f ish3.pgm"). 

The noise was added "artificially" by a program which scan every bit in 
the image and inverts it with a given probability. The probabilies used were 
10% and 15%. The modified files were named fish. lO.pgm, f ish2. lO.pgm, 
f ish3 . 10 .pgm and fish. 15 .pgm, f ish2 . 15 .pgm, f ish3 . 15 .pgm, respectively. 

Each image was restored in all bit planes using 3 values of (3, namely 0.1, 
0.3, 0.5. An aditional restoration, just on the most significant bit plane were 
made using the estimated (3 as indicated in the appendix. The values are 
presented in the tables ||, ^| and [|. 

We did not try to define a metric to measure the quality of the recon- 
struction, relying on a subjective analysis. Besides the removal of noise, it 
was observed an improvement in the shades and (consequently) in the third 
dimension perception, even if some blurring is introduced. 

These and other examples are available directly from the authors or at 
http : / /www . ime . usp . br/~gubi/MAP/index . html 



File 


Beta 

0.10 0.30 0.50 


fish. lO.pgm 
fish. 15. pgm 
f ish2. lO.pgm 
f ish2. 15. pgm 
fish3. lO.pgm 
f ish3. 15. pgm 


6.40 17.77 16.78 
6.05 14.38 19.54 
2.81 4.45 6.35 
2.65 4.38 5.86 
1.38 2.43 3.20 
1.34 2.37 3.02 



Table 1: Restoration running times for three sample images with two noise 
levels (10% and 15% per bit) 
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Table 2: Sample image (fish.pbm) 
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3 MAP estimator and networks 



In this section we present some capacitated network ideas which are used to 
solve the two-color maximization problem described before as discovered by 
Greg, Porteous and Seheult [GPS] and use them to prove Proposition. 

As mentioned before with the choice of H given by (fP2| ) the maximization 
in (j^) is decomposed into K two-color problems. Therefore we assume in this 
section that we are in the binary color case K = 1. 

The MAP estimator given Y is the image X e {0, 1} A which gives the 
maximum of 

L a (X\Y) = aJ21 {Xi=Yi} (X) + £ l {Xi =x j} (X) (15) 

ieA <i,j> 

where a = h/ (5. 

A network iV is a graph G = (V, E) , where V is a finite set of vertices 
and E is a set of couples of vertices, with a capacity c(e) > associated with 
each edge e e E. 

We define networks on the set of vertices V, given by the sites in A plus 
two extra ones denoted by s (source) and t (sink), that is 

V = AU {s} U {*}. 

For the set of arcs we take 

E = (Ui{(s,i)}) U (Ui{(i,t)}) U (U <id> {(i,j)}) 

where the first two unions are taken on i G A and the last one is taken over 
pairs of nearest neighbor sites in A. 

Given the observed image Y_ and a real number a we define the capacities 
of the network N a (Y_) as follows. 

If Y(i) = 1 we set cy, a (s,i) = a as its capacity, otherwise set cy )0l (i, t) = 
a; to each arc e = of neighbor sites in A we associate cy, a (e) = 1. All 
other arcs have capacity zero. 

For each image X_ let 

A(X) = {s} U {% e A : X{%) = 1} 
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B(X) = {t} u {i E A : X(i) = 0}. 
These two sets define a cut of the network 

CQQ = eE-.ie A(X),j e B(X)}. 

Notice that the cut C(X) consists of a set of arcs whose removal (cut) makes 
it impossible to find a path going from s to t through arcs with non-zero 
capacity. If we now define the capacity of the cut C(2Q by the quantity 

CV>(X)= £ cvAhj) (16) 
(ij)ec(x) 

it is very simple to check that 

L a (X\Y)=a-Cy, a (X) } 

where a is a constant which does not depend on X_. Therefore to find the 
MAP estimator one has to find the cut which minimizes flip]): the so called 
minimum cut [FF]. 

Ford and Fulkerson showed that the value of the capacity of the minimum 
cut is equal to the maximum flow through the network from source to sink. 
Recall that a flow f in a network N (on G(V, E) with capacity {c(e)} e€ £) from 
s to t is a collection of real numbers {f(e)} ee £, where f(e) can be thought 
as the amount of fluid per unit time going through the pipeline e, such that 
the flow on each arch does not exceed its capacity, < f(e) < c(e), for all 
e £ E, and such that the flow is conserved 

£f(i,j)-£f(.7,i) = o (17) 

jev jev 

for all i s. 

Well known fast algorithms to find this maximum flow exist and therefore 
the problem is solved in practice. The algorithm presented in this paper uses 
this method to reconstruct each one of the (8 or 24) layers. 

Proof of the Proposition. To simplify the notation write X_ = E a (Y_). Let 
a' > a. The networks defined by the pair (Y_, a), denoted by N a (Y_), and by 
( X , a'), N a i(X_), can differ only in the capacities assigned to each arc. 
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The proposition asserts that the cut defined by X_ also has minimum 
capacity in N a i(X_), that is we want to prove that 

C^(X)<C Laf (Z), (18) 

for any image Z G 

Fix Z and let Ej = C(X)\ C(Z), E u — C(X_) H C(Z) and E IU = C(Z)\ 
C (X) . By definition of MAP estimator we have 

Cy, a (l) < Cy >a (Z), (19) 

which implies 

< DyAEiii). (20) 

where, for a set of edges i<7, 

Dy,a( E ) = H Cy iQ (e) 

Therefore to check (|18|) it is enough to verify 

D^JEj) < Dy, a (Ei) (21) 

D Lal {E UI )) > Dy^Ejjj). (22) 

We start with inequality fl2"T|). Assume e G Ej. If e is an external edge 
(connecting some site i with either the source s or the sink t), then, since 
e G C(A), it must be c^ a ,(e) = which is less or equal than cy, a (e). On 
the other hand, if e is an internal edge (i.e. e = for i and j nearest 

neighbors in A) then its capacity equals 1 in both cases. This concludes the 
proof of (gT|). 

Suppose now e G E HI . Inequality fl2"2|) follows from the observation that 

C X a'( e ) CaI1 n °t ^ e zero - □ 

4 Final Remarks 

The different approaches used in image reconstruction are based in quite dif- 
ferent set of theoretical ideas and it is not clear how to compare their results. 
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One possible measure for the quality of the reconstruction, used in [FP] 
to compare 9 algorithms, is to evaluate the proportion of pixels classified 
correctly. Since our main goal here was to present a method for the recon- 
struction of multicolor images we leave the comparison with other methods 
for future work. In any case, since we are working with exact MAP's for the 
chosen Hamiltonian, our method will be as good (and as bad) as the usual 
two colors MAP estimators, regarding the proportion of bits reconstructed. 

For the usual MAP reconstruction problem in the multicolor case no 
fast algorithm is known ([FFG]). For a probabilistic approach via Simulated 
Annealing in order to get the exact estimator one needs to decrease very 
slowly some parameter while the computation goes on and thus needs a 
prohibitively large amount of time [Gi]. 

One possibility is to accept approximated MAP estimators which can be 
obtained fast enough. One can do this with simulated annealing by updating 
the parameter fast enough but then we lose control on how close the approx- 
imation is to the exact one. An approximate of the MAP estimator with a 
probabilistic analysis of the error in the three color case was developed in 
[FFG]. Another approach can be found in [J]. 

Our approach is not completely Bayesian as we also consider a situation 
on which the a priori measure has some unknown parameters but do not 
assume any prior knowledge about them. In the Appendix we describe a 
method to estimate these parameters from the picture itself using classical 
statistical methods, as proposed by Frigessi and Piccioni [FP]. In those cases 
we verify that plugging those estimated values into the formulas used to get 
the MAP estimator provides a reconstruction which appears to be the best. 
As mentioned before we do not try to quantify this. 

Ricardo Maronna proposed that instead of estimating (3 and h one could 
look for the a = (3/h that maximizes (0). The estimator for the true image 
will then be the X_ which realizes this maximum with the best a. It is 
not clear yet how to justify this theoretically. Another alternative would 
be to choose the uniform distribution on an appropriate range as a priori 
distribution for e and f3. However the computation of the MAP estimator in 
this case seems prohibitive. 

As an experimental observation we remark that reconstructing only the 
first layer of a dirty image (and leaving the others as they are) gives a quite 
good visual result. A possible explanation of this fact is that, as a conse- 
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quence of the binary decomposition, each layer is "half" as important as the 
previous one. 

Looking for algorithms which give smooth solutions we propose the fol- 
lowing hierarchical procedure. First consider layer 1 and find Z 1 , as before 
for X 1 . Then given layer Z e , for £ — 1, . . . , k — 1, define Z k as the (binary) 
image that maximizes 

k 

L a {X k \Y) = a ^l {xhY H } UO + E II l { xf=xf}00- (23) 

iGA <i,j> £=1 

In words, this algorithms will try to get the same value for two neighbors 
in layer k if these neighbors have the same value for all previous layers. If 
not they are not coupled. Visual realizations of this algorithm give also good 
results. Although there is no Hamiltonian for this model, the algorithm is 
well defined. Moreover each layer corresponds to the so called diluted Ising 
model. 

Assume that we get N independent samples of the same image Hi, ... , Y_ N . 
This means that the original realization X_ of the image is the same but the 
noises are independent. A generalization of our approach deals with this 
problem in the following way. We modify the capacities associated to the 
graph G. To each arc (s, i) connected to the source we associate the capacity 

and then to each arc (i, t) connected to the sink we associate the capacity 

"(EW.i- 1 )) • 

In the same vein one could use the remaining two colors to get information 
about the color being reconstructed. 

5 Appendix: Estimation of /3 and e 

In this section we apply ideas from Frigessi and Piccioni ([FP]) which exploit 
well known (but highly nontrivial) results from statistical mechanics ([R], 
[MW]) to obtain estimators for the parameters f3 and e = (1 + e fe ) _1 . 
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The result is a program called estima that has as input a multicolor 
image and returns one estimated (3 and one estimated e for each layer. 

The measure defined in (§) with H given by @ in the two-color case 
(Ising model) in the thermodynamic limit gives rise to a translation invariant 
measure which may be ergodic or not according to the value of (3. 

This choice of H does not favor zeros or ones. This symmetry is perhaps 
easier to see if we represent a configuration as an element S_ in {— 1, 
with —1 replacing zeros. This is the usual Ising notation while the one using 
zeros and ones is known in statistical mechanics as lattice gas representation. 
In the situation considered here they are equivalent. 

More precisely in if F G {0, 1} A let S_ be the corresponding Ising config- 
uration given by 

S={S i ,ieA:S i = 2Y i -l} (24) 
And instead of Hamiltonian (Bf) we use 

H P (S) = - £ SiSj (25) 

<i,j> 

The corresponding infinite volume limit measure is not ergodic if (3 is 
smaller than some (known) critical value (3 C . In this case the limit measure is 
a mixture with equal weights of two measures //^ and [ip, the first favoring 
configurations with more +l's than — l's and the other favoring configura- 
tions with more —l's than +l's. If (3 > (3 C then fit = If we denote by 
Et(f) ( E7(f) ) the expected value of a function / defined on { — 1,+1}^ 2 
the symmetry between \x~p and fi^ imply 

EfiSiSj) = E^S^) = r„(\i-j\) (26) 

r p{\i ~~ j\) is called the two point correlation function for the infinite 
volume system at inverse temperature (3. We will need two of those correla- 
tion functions in what follows: ^(1), for nearest neighbors, and rg(v^2) for 
neighbors along the diagonal on the lattice. 

Suppose S_ corresponds (in the Ising notation) to the original image and 
R corresponds to the observed image after noise. Under the hypothesis on 
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the noise 



E+(R i R j )=E+(S i S j )(l-2e) 2 (27) 

Therefore the ratio between 773(1) and 773 (\/2) computed for the observed 
image depends only on /3. Call this ratio 4>(/3). 

If Re { — 1, 1} A corresponds to an observed two-color picture [FP] found 
a sequence (/3, e) of consistent estimators for (f3, e) given by 

P = <P 

and 
with 

G'(R) = n ^ R ' (30) 

where for a =1 the sum is taken over pairs of nearest neighbor sites along 
the lattice directions in A , the interior of A, and for a = 2 the sum is over 
neighbor sites along the two lattice diagonals again in the interior of A. 

Expressions for 773(1) and rp(y/2) are complicated involving elliptic in- 
tegrals with different formulas for (3 smaller, equal and larger than j3 c and 
the function _1 must be computed numerically. Our program estima finds 
these estimators from a observed image. 
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